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NOMENCLATURE 


S^mbo 1 


■Jd 


excitation  amplitude 
post-yield  stiffness 
energy  of  dissipation 
a  function  of  random  variables 
permanent  set 

stiffness  in  the  elastic  range 

system  mass 

maximun  displacement, 

tilOA 

restoring  force 
time 

the  time  the  system  reaches  its  maximun  displacement, 

iTiqLX 

the  time  the  system  first  reaches  its  yield  displacement,  u 7 

yield  displacement 

response  displacement 

maximum  displacement  response,  R 

excitation  decay  rate 

damping  ratio 

mean 

random  variable 
correlation  coefficient 
standard  deviation 
random  variable 

natural  frequency  (small  displacement) 
damped  frequency  (small  displacement) 
post-yield  frequency  (■  /A 7/m) 


v 


I 

\ 

1.0  ^JNTRCOUCTION 

An  above  ground  explosion  generates  a  shock  wave  in  air,  or  airblast, 
and  is  accompanied  by  some  duration  of  strong  wind.  Especially,  the  air- 
blast  from  a  nuclear  detonation  can  cause  extremely  high  air  pressures  and 
has  a  relatively  long  period  of  duration.  When  structures  under  design  may 
be  subjected  to  this  sort  of  loading  condition,  it  is  necessary  to  analyze 
the  behaviors  of  the  stuctures  when  subjected  to  this  kind  of  input  to 
determine  whether  the  structures  can  survive  or  not. — ^  p  I^IS 

When  an  airblast  produces  a  high  air  pressure  and  has  a  long  period  of 
duration.  It  will  cause  many  structures  to  have  an  extreme  response.  , If 
the  structure  can  execute  a  response  that  has  plastic  deformation,  the 
designer  must  have  some  technique  to  analyze  Inelastic  response  in  the 
design  process.  The  designer  needs  the  capability  to  predict  such  measures 
of  structural  response  as  the  plastic  deformation  and  maximum  displacement. 

When  a  structure  executes  an  extreme  response,  it  may  behave  inelas- 
tlcally,  and  during  the  response,  offsets  from  the  initial  elastic  configu¬ 
ration  may  occur..  When  these  offsets  are  large,  residual  deformation  in 
the  structure  may  exist  at  the  completion  of  structural  motion,  and  the 
magnitudes  of  these  residual  deformations  may  be  important.  When  they  are 
Important,  It  is  desirable  to  Include  the  potential  for  these  residual 
deformations  In  the  mathematical  model  of  the  structure  so  that  they  may  be 
predicted  In  the  structural  analysis.  The  models  used  In  References,  [1], 
C2],  [3]  do  account  for  some  features  of  Inelastic  response,  but  do  not  , 
permit  permanent  offset  (residual  deformation).  The  model  used  in  the 
current  investigation  does  account  for  permanent  offset. 

The  properties  of  nuclear  and  high  explosive  airb lasts  are  random. 

The  reason  Is  that  the  material  and  geometric  properties  of  blast  sources 
are  random.  In  view  of  this,  an  airblast  measurement  is  never  repeated' 
even  when  nominally  identical  blast  sources  are  used.  An  airblast  signal 
usually  consists  of  a  rapidly  Increasing  pressure  followed  by  a  gradual 
decay.  Two  properties  of  the,  pressure  wave  are  of  great  concern.  The 
first  is  the  amplitude  of  the  airblast;  the  second  is.  the  decay  rate.  The 
first  property  describes  the  highest  air  pressure  the  airblast  produces 
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throughout  the  time  history.  In  general,  the. rise  time  is  so  brief  that 
the  peak  amplitude  Is  assumed  to  occur  at  the  beginning  of  explosion.  The 
second  property  deals  with  the  rate  at  which  the  airblast  dies  out.  In 
this  report,  these  two  properties  are  considered  random. 

In  most  applications  the  structure  itself  should  also  be  considered 
random.  The  material  properties  of  the  structure  determine  the  elastic 
stiffness,  yield  stiffness,  yield  displacement,  etc.  Conseouently,  the 
natural  frequency  Is  affected  strongly  by  the  material  properties.  When 
the  material  of  the  structure  is  considered  to  be  random,  then  the  proper¬ 
ties  of  the  structure  must  also  be  considered  as  random.  In  this  report, 
the  material  properties  considered  as  random  are  natural  frequency,  yield 
stiffness,  yield  displacement  and  damping  ratio.  These  four  properties  are 
considered  random  because  they  affect  the  response  significantly. 

The  previous  discussion  points  out  that  both  the  excitation  and  struc¬ 
tural  system  can  be  regarded  as  random;  in  this  study  both  sources  of 
randomness  are  considered.  The  response  of  the  system,  therefore,  can  be 
regarded  as  a  response  random  process,  and  there  is  a  relationship  that 
characterizes  the  response  In  terms  of  the  excl talon  and  system  random 
variables. 

These  random  variables  Include  Input  random  variables  which  are  peak 
amplitude  of  air  pressure  and  decay  rate,  and  system  random  variables 
which  are  natural  frequency,  yield  stiffness,  yield  displacement,  and  damp¬ 
ing  ratio.  It  Is  desirable  to  develop  a  method  to  establish  the  relation¬ 
ship  among  the  Ifiput  characteristics,  the  system  characteristics,  and  the 
response  characteristics.  The  response  random  process  Is  then  character¬ 
ized  In  terms  of  the  relationship. 

Many  measures  of  structural  response  concern  the  designer.  Some 
measures  of  special  importance  are  plastic  deformation,  maximum  displace¬ 
ment  response,  and  energy  dissipated  by  the  system  spring,,  or  by  the 
damper.  The  plastic  deformation,  which  Is  the  offset  from  the  static 
configuration,  must  be  controlled  in  structures.  The  maximum  displacement 
response  must  be  less  than  the  design  value.  The  energy  dissipated  must  v 
less  than  an  allowable  value  since  It  Is  related  to  the  fatigue  and  failur- 
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of  a  structure.  When  the  probabilistic  character  of  the  input  excitaion 
and  system  random  variables  are  known,  and  when  the  relationship  between 
measures  of  response  and  the  input  and  system  variables  has  been  estab¬ 
lished,  then  it  is  possible  to  specify  the  response  random  process. 

Certain  probabi listic  measures  of  structural  response  are  obtained.  Spe¬ 
cifically,  the  mean  and  variance  of  permanent  offset,  peak  response,  and 
energy  dis sip ted  in  an  Inelastic  structure  are  obtined.  The  results  can  be 
used  to  determine  whether  design  criteria  are  satisfied. 

In  past  investigations,  several  mathematical  models  have  been  used 
with  structural  identification  procedures  to  study  structural  behavior. 
Lutes  and  Hseih  [1]  used  linear,  third  and  fourth  order  equations  to 
approximately  simulate  the  response  of  an  SDr  hysteretic  system.  They 
found  that  the  spectral  density,  or  the  autocorrelation  functon  of  the 
linear  and  hysteretic  systems  can  be  matched.  Toussi  and  Yao  [2]  used  the 
response  of  a  ten-story  reinforced  concrete  test  structure  as  a  means  to 
study  the  feasibility  and  practi callity  of  a  hysteresis  identification 
technique.  They  considered  the  hysteresis  to  be  a  measure  of  damage.  Yao, 
Toussi  and  Sozen  presented  the  state-of-the  art  on  damage  assessment  of 
existing  structures  in  Reference  [3].  Morrison  and  Paez  [4]  proposed  a 
procedure  for  predicting  the  probability  of  survival  of  inelastic  struc¬ 
tures  excited  by  blast.  In  this  study  the  failure  criterion  was  related  to 
peak  strain  in  the  structural  repsonse.  A  method  for  establishing  the 
probabilistic  character  of  nonlinear  structural  repsonse  was  developed. 
Rohanl  [5]  presented  a  probabilistic  solution  of  one-dimensional  wave 
propagation  phenomenon  In  earth  media  by  using  a  deterministic  model. 
3ennett  and  Paez  [53  provided  a  means  for  relating  uncertainties  of  Input 
parameters  to  the  uncertainty  in  response  in  a  reinforced  concrete  slab. 

The  method  used  In  the  last  three  references  Is  adopted  In'  this  investiga¬ 
tion  to  characterize  nonlinear  .measures  of  inelastic  structure  response. 

The  objectives  of  this  Investigation  are,  first,  to  develop  a  model  to 
define  the  permanent  offset  In  an  inelastic  structure.  Next,  we  develop 
the  relationship  among  input  characteristic,  system  characteristic  and 
response  characteristic  through  which  the  response  random  process  Is 
established.  Based  on  this  response  random  process,  and  given  the  moments 


of  input  and  system  variables,  the  response  characteristics  are  computed. 
This  provides  a  means  to  analyze  the  structural-  response  when  the  structure 
is  subjected  to  blast  type  loading  and  ooth  the  input  and  the  system  are 
considered  random. 

The  structure  considered  here  is  a  single-degree-of -freedom  (SDF) 
bilinear  hysteretic  system  which  possesses  an  elastic  stiffness  and  yield 
stiffness.  The  natural  frequency  and  yield  natural  frequency  are  therefore 
defined.  The  hysteresis  loop  of  a  spring,  strongly  influences  the  energy 
dissipated  in  a  structure;  this  will  be  discussed  in  detail.  The  airblast 
input  is  modeled  by  a  decaying  exponential  function.  The  value  of  peak  air 
pressure  and  decay  rata  of  the  input  are  obtained  from  che  A.ir  Force  Design 
Manual  [7]. 


2.0  THE  ANALYSIS  OF  MAXIMUM  RESPONSE 

Tne  maximum  response  of  a  structure  wnen  subjected  cc  blast  type  load 
can  be  predicted  through  basic  structure  response  analysis.  Intuitively, 
the  maximum  response  must  relate  to  the  input  parameters,  which  are  peak 
air  pressures  and  decay  rate,  and  system  parameters,  which  are  natural 
frequency,  damping  ratio,  etc.  All  these  properties  affect  the  response. 
Specifically,  we  can  establish  a  mathematical  function  which  governs  the 
maximum  '•eponse.  The  response  random  process  can  be  simply  expressed. 

Given  the  statistical  moments  of  the  input  and  system  parameters,  the 
statistical  moments  of  the  responses  can  be  computed.  The  response 
considered  here  can  be  inelastic. 

It  is  assumed  that  the  systan  under  consideration  has  significant 
elastic  and  yield  stiffness  so  that  the  bilinear  hysteretic  system  can  be 
used  as  the  structural  model.  The  offset  function  for  blast  type  input 
wil1  behave  as  in  Figure  2-1. 

Again  consider  an  SDF,  bilinear  hysteretic  system,  with  small  damping 
ratio,  c,  system  mass,  m,  elastic  stiffness,  k.  The  equation  of  motion  is 
given  by 

Z  *  2 ;unZ  +  ^(z)  »  A'at  t  >  Q  (?-l) 

Assume  that  the  system  starts  from  rest,  z(0)  *  0  and  z'(0)  ■  0.  Let  t7 
denote  the  time  the  system  first  reaches  its  yield  displacement,  U7. 

Within  time  t  <  t^,  the.  system  remains  linear,  so  that  Equation  (2-1) 
can  be  written  as 

z  f  2c  unz  +  mjjz  *  Ae’3*’  for  t  <  ty  (2-2) 

Because  the  system  remains  linear  during  the  initial  portion  of 
response,  the  above  equation  can  be  solved  Immediately  by  using  convolution 
Integral 

4» 

z(t)  •  /  h(t)  X  (t  -  ?)  dr  (2-3) 

0 


where  h(t)  Is  the  linear  system  impulse  response  function  and  x(t)  is  the 
input  forcing  function.  The  impulse  response  function  is 


h(t) 


.  sin  Udt) 


(2-4) 


where 

z(t) 


“d  *  “n  ^T-d  x(t)  «  Ae“at.  Carry  out  the  integral;  then 


A  ( 

.** - 2  1 . 21  )e  C(«  *  sin(o>dt)  -  w^ccsioj^t)]  + 

(«  *  5^)  +  “dj* 

fdr  0  <  t  £  tj 


(2-5) 


At  time  t"  »  t7,  the  displacement  response  should  be  equal  to  u7. 


z(t7) 


-Ctu„t 


7  a>d  [(a  -  cu»n)^  +  o»jj]l 


e  n  n[(a  -  Ca^isinfcj^t^) 


COS(a»dty)  + 


•=t7) 

»de  | 


(2-6) 


The  time  ty  at  'which  yielding  first  occurs  must  be  determined  by 
solving  the  above  equation  for  ty.  Since  q  and  t7  are  both  small; 

.•“‘"'7  1. 


*  p;  o>-,[P2  +  «?] 


Simplify  the  notation  by  letting  a  -  -  h»  a>du*  •  «.d. 

«  x^;  A/x^  *  Xg.  Rewrite  Equation  (2-6)  using  the  above  notation. 


-  u7  (2-7) 

The  value  of  «  is  about  from  0.5  -  6.0.  This*  is  usually  much  smaller 
than  the  value  of  <jn;  hence, 

p  •  sin(udt7)  «  udcos(udt7) 
and  the  quantity  on  the  left  is  negligible. 


<2  .  p  Sin(a>dt7)  -  ^  COS(adt7)  +  ad 


•at- 


Using  Taylor's  expansion  to  express  e  and  cos'Cu^t^)  and  neglect¬ 
ing  the  high  order  term,  t-/  can  finally  be  evaluated 

t?  »  1  (B  +^3  +  4C  where  3  »  2o/(<!|  +  a2),  C  *  2(u?)/A.  (2-8) 

After  solving -for  t7,  the  velocity  at  t  »  tj  can  be  computed  immedi- 
•  • 

ately.  Let  z(t7)  *  z^. 

To  this  point,  the  only  assumption  used  in  determining  the  time  of 
first  yield  and  the  velocity  at  that  time  is  that  the  system  viscous 
damping  factor  is  relatively  small.  No  assumptions  regarding  the  input 
level  or  yield  level  have  been  used. 

Next,  we  consider  the  phase  of  the  response  when  yielding  starts  to 
occur.  Let  t,,,^  be  the  time  the  response  reaches  its  maximum 


displacement  value.  Ouring  the  time  interval  (t7,  t^a*)  the  system  has 
entered  the  plastic  region  and  permanent  offset  starts  to  occur.  If  -we  let 
A7  denote  the  yield  stiffness,  then  for  time  <  tj,  the  restoring 

force  function  can  be  written  as  (refer  to  Fig.  2-2): 

I  ♦  Z  ♦  ; 

|CA 7Z  +  (k  -  Ay)  u?]  ■  Ae  0  , 

tmaY  >  t  >  t 
max  —  — 

(2-9) 

or 

2  +  2«anZ  + 

;z  »  Ae"at  -  (ujj  -  «y)u7  , 

w  >  * >c 

(2-10). 

2 

where  a*  *  A^/m 

. 

The  solution  i 

'or  z(t)  can  be  expressed 

"ta-t 

Z(t)  -  e  n  ( 

c ^cos(u^t)  *  cz  sin  ;dt)  +  zp  , 

s  i ‘  <  W 

(2-11) 

where  •  •J(ky'm) 

*  (^n)2 

. 

fe  ^TCT^T-'VW? W.W 


2p  *  ql  +  q2  e 

qx  *  [(ky/m  -  zy  m/ky 

^2  *  A/[a^  -  2qaina  +  (ky^m)] 

By  using  q^  and  qj  in  Equation  2-11,  the  constants  and  Cg  can  be 

solved  by  using  as  initial  condition,  z(tj)  *  u^  and  z(t^)  »  z-j  which  were 
obtained  in  the  previous  discussion. 

Since  the  displacement  response  has  been  determined  in  Equation 
(2-11),  the  maximum  value  can  be  obtained  by  maximizing  Equation  (2-11). 
The  maximum  occurs  when 


z  <‘m»>  ’  0 


(2-12) 


The  ^ax  can  be  obtained  from  Equation  (2-12).  The  corresponding 
then  can  be  evaluated  by  using  t^g*  in  Equation  (2-11). 


10 


3.0  PROBABILISTIC  ANALYSIS  OF  RESPONSE 

In  the  previous  section,  the  relationship  among  the  maximun  response 
and  input  parameters  and  system  parameters  was  established.  The  response 
random  process  is  then  characterized  in  terms  of  the  established  relation¬ 
ship.  In  many  cases,  the  established  function  may  be  complicated.  In 
order  to  analyze  some  measures  of  the  response,  a  method  using  Taylor's 
expansion  Is  available  and  has  been  widely  used.  A  brief  discussion  of 
that  method  is  presented  here. 

Let  R  denote  a  random  variable  that  is  a  measure  of  the  structural 
response  random  process.  It  is  assuned  that  R  is  a  function  of  n  random 
‘  variables.  These  n  random  variables  are  the  parameters  of  the  input  and 
the  structural  system.  Let  the  functional  relation  be  denoted 


R  *  9(yi»  Y2» . Yn) 


(3.1) 


The  Y-j*  i*i...n  are  the  random  variable  with  mean  values  ,  i*l,...n. 

The  function  g  can  be  expanded  in  a  Taylor  series  about  the  means  of  the 
parameters: 


R  *  9(u  »  u  . u  )  +  2Z  (y-{  -  u*  )  + 

where  }  is  the  vector  of  mean  values  of  the  y^. 


(3-2) 


The  moments  of  R  can  be  evaluated  by  using  sane  portion  of  the  series 
in  Equation  (3-2).  For  example,  the  series  can  be  truncated  following  the 
linear  terms.  When  this  is  done,  the  mean  of  R  is 


SW  *  9<V  V  ••••-,  ) 


(3-3) 


This  mathematic  a]  expression  shows  that  the  mean  value  of  the  depen¬ 
dent  random  variable,  R,  is  closely  related  to  the  mean  values  of  the 
underlying  random  variables,  >  and  the  function  expression  g(u^  ,  ^  .. 

K  )*  1 

Tn 

The  variance  of  R  can  be  computed  as: 


Var[R] 


E[R2]  -  E[R]2 


*  E J?  (uyp  ,JY2’  *  2J29(“y1’uy2’  •••uYn)  • 

4  ] 

)%r 


n  n 

£  2  pij 

1*1  j-i  1  1  3Yi 


W 


(3-4) 


The  a.  and  ^  represent  the  standard  deviations  of  y^  30(1  Tj  and  here 

represents  the  coefficient  of  correlation  between  y^  and  y,. 

•+ 

The  analysis  technique  outlined  above  can  be  used  in  a  wide  variety  of 
cases;  especially,  it  is  useful  in  the  characteristics  of  the  response  of  a 
nonlinear  system.  When  the  moments  of  the  underlying  random  variables  are 
known  (means,  variances,  and  covariances) ,  then  it  is  easy  to  apply.  But 
there  are  some  cases  where  tne  random  variable  R  depends  on  some  underlying 
random  variables  whose  characteristics  are  unknown,  a  priori,  as  well  as 
random  variables  whose  characteristics  are  known.  For  example,  we  may  be 
interested  in  establishing  the  moments  of  energy  dissipated  during 
inelastical  response.  This  quantity  can  be  expressed  as  a  function  of 


•r»«- 


•X 


*vr  r  *v*.*v  cv  t-v  tr*  cttittc 


input  parameters,  system  parameters  and  the  peak  response.  The  moments  of 
the  input  and  system  parameters  are  probably  known,  a  priori,  but  the 
moments  of  peak  response  certainly  are  not.  To  account  for  this,  the  first 
phase  of  an  investigation  can  be  aimed  at  analyzing  peak  response  and  its 
statistical  relation  to  the  input  and  system  parameters.  After  this  is 
done,  the  energy  dissipated  can  be  analyzed  using  the  general  approach 
outlined  above.  In  general,  any  sequence  of  analyses  like  this  one  can  be 
executed.  At  each  step,  a  new  random  variable  characterizing  the  response 
can  be  introduced. 

For  example,  the  analysis  referred  to  above  can  be  accomplished  in  the 
following  way.  Let  R  be  a  function  of  n  random  variables,  si  1*1,2,... m, 
m**l,m+2,....n.  The  moments  and  cross  moments  of  the  ai,  1*1, 2, ....m,  are 
known,  a  priori,  and  the  moments  and  cross  moments  of  the  f?i,i*m+l,....n, 
are  not  known.  '  However,  the  moments  of  the  Si,iam+l,....n,  can  be  deter¬ 
mined  from  the  8i,i*l,....ffl.  Let  Si,i*m+l,....n  be  expressed  as 
follows: 

31*  h^s^Sg . Sm),  1*m+l,...n  (3-5) 

The  Taylor  expansion  for  the  above  expression  can  be  writtenr  as  in  Equa¬ 
tion  (3-2).  The  series  can  be  truncated  at  some  point  and  the  result  can 
be  used  to  approximate  the  moments  and  cross  moments  for  all  the 
8^  ,1*m+l, . . .n.  For  example,  when  the  series  are  truncated  following  the 
linear  term,  the  covariance  between  two  of  the  s^'s,  one  from  the  group 

1*1,... m,  and  the  other  from  the  group,  i*m+l,...n,  can  be  expressed  as 


Cov(8pSq) 


The  means  and  variances  of  the  sj,  1»m+l,....n  are  obtained  in  the 
same  way  as  Equations  (3-3)  and  (3-4).  Given  these  moments  and  cross 
J  moments,  the  mean  and  variance  of  R  can  be  obtained  using  Equation  (3-3) 

;  and  (3-4),  as  they  are  vcitten. 

;  The  response  randan  process  can  be  accurately  characterized  using 

.three  measures:  the  permanent  offset,,  the  max Imur  response,  and  the  energy 
|  dissipated.  Using  the  techniques  described  above,  the  means  and  variances 

of  each  random  variable  can  be  computed.  The  following  section  will 
present  more  detail  about  those  randan  variables. 

I  3.1  The  Maximum  Response  Random  Process 

From  Section  2.0,  the  maxlmun  displacement  response  is  governed  by  5 
random  parameters  and  can  be  represented  as: 

l  3  •  9(rj»  Y2 — ^6) 
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where  R  *  maximum  displacement  response 

*  input  peak  air  pressures  (A),  y2  3  decay  rate  (a) 

Y3  *  natural  frequency  Un),  Y4  “-^Ay/m 

Yg  *  yield  displacement  u7,  Yg  a  damping  ratio  ($) 

Using  Equations  (3-3),  (3-4),  the  mean  and  variance  of  R  can  then  be 
computed. 

3.2  Tha  Permanent  Offset  Random  Process 

The  relation  between  permanent  offset,  H,  and  maximum  displacement,  R, 
i  s : 

A? 

H  -  (R  -  uyjd  -/)  (3-8) 

So  we  can  let 

H  «  h(3x,  s2,  83.  84) 

•where  3^  *  un,  s2  *  S3  ■  U7,  84  *  R:  Then  the  mean  and  variance 

of  H  can  be  computed  using  the  technique  as  stated  above. 

3.3  The  Energy  Dissipated  Random  Process 

The  energy  dissipated  by  a  bilinear  hystaretlc  system  Includes  two 
parts.  The  first  part  Is  the  energy  dissipated  by  the  spring.  The  second 
part  Is  the  energy  dissipated  by  the  damper.  The  system  will  reach  Its 
maximum  offset  within  the  first  cycle.  Then  the  energy  Is  dissipated  by 
the  spring  only  during  this  time.  After  the  .system  reaches  Its  maximum 
offset,  then  the  system  remains  linear  in  a  new  equilibrium  configuration. 
The  energy  is  then  dissipated  by  the  damper  only.  The  energy  dissipated  by 
a  linear  system  Is  related  to  the  damper  only  and  has  already  been  studied 
(Reference  8).  In  this  report,  only  energy  dissipated  by  the  spring  is 
considered. 

Refer  to  Figure  (2-2).  The  energy  dissipated  by  the  spring  can  be 
computed  as: 
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Ed  *  *5  *  k  *  u7  +  -5  (2aax  -  u7)A?  +  k  .  u7(2max  -  «7)  -7  k(2max 

therefore  £d  *  ^u7  +  7  ’  u7)2  +  “nu7(2max  *  u7} 

*  T^max  "  u7)“y/“n  +  VJ* 

In  other  words,  the  energy  dissipated  can  be  expressed  as: 


-H)2 

(3-9) 


Ed  *  s(ni,  n2t  n3*  nij.  (3-10) 

■where  we  denote  m  *  Wj,,  n*  *  ay  nj  *  uy,  n*  *  z^.  Again  using 

Equation  (3-3),  (3-4),  and  (3-6),  the  mean  and  variance  of  energy 
dissipated  can  be  computed. 


4.0  NUMERICAL  EXAMPLE' 

In  Section  2.0  an  approach  for  computing  the  maximum  displacement  for 
a  bilinear  hysteretic  system  subjected  to  blast  type  load  was  developed. 

In  Section  3.0  an  approach  for  the  probabilistic  analysis  of  the  response 
was  developed.  A  computer  program,  called  YARE.F  Is  used  in  this  section 
to  calculate  the  means  and  variances  of  the  response,  which  includes 
maximum  displacement,  permanent  offset  and  energy  dissipated  by  the  sprin. 

In  Section  3.0  the  technique  used  to  evaluate  the  variance  of  a 
response  random  variable,  which  depends  on  n  underlying  random  variables 
and  whose  statistical  moments  are  known,  refers  to  the  partial  derivatives 
with  respect  to  the  underlying  random  variables.  In  order  to  calculate  the 
partial  derivatives  with  respect  to  the  underlying  variables,  an  ' 
approximated  method  Is  used  here,  let  R  be  a  random  variable  which  depends 
on  n  underlying  variables  y^,  1*1, ....n  (for  example,  maximum 
displacement  which  Is  expressed  in  Equation  3-1.).  Then 

R  ■  9(y1*  Y2»»»»Yn)  (4-1) 

and  the  partial  derivative  of  R  with  respoect  to  y-f  can  be  approximated 
12.  .  g(rl’  y2  Y1  +  ^1  •••V’9  {y1»  r2  "‘T1  *  ATi  —  rn}(1  2) 

This  approximation  can  be  computed  even  though  g(r^»Y2»  •  ••  Yn)  may 

not  be  a  function  that  can  be  written  explicitly.  Writing  g(y^»Y2  •••Yn) 

here  Implies  that  a  computer  program  can  be  run  using  the  y^,  1*1,... n,  as 

Inputs  to  obtain  a  result  R.  The  program  can  be  run  twice  to  obtain  the 
above  approximation,  using  aoove  approximation  In  the  program  YARE.F  the 
first  order  partial  derivative  with  respect  to  any  given  underlying  random 
variable  can  be  calculated.  The  maximum  displacement,  permanent  offset  and 
energy  dissipated  are  considered  in  the  following  section. 

4.1  The  Maximum  Displacement 

It  Is  mentioned  In  Section  3.1  that  the  maxlmun  displacement  is 
rel»ted  to  six  random  parameters.  The  functional  expression  is 
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R  *  9(n»  T2*  Y3>  Yu.  Yu.  Y5 .  Yg) 

where  R  is  maximun  displacement,  yi  -  input  peak  air  pressure 

Y2  *  decay  rate  of  input  force,  Y3  3  natural  frequency,  u>n 

Yu  *  /A77m  Y5  3  damping  ratio, 

Y6  3  yield  displacement,  u7 

The  means  and  standard  deviations  of  the  above  six  randait  variables 
are  given  in  Table  4.1. 

Table  4.1 

The  Means  and  Standard  Deviation  of  Yj,  i*l,6 

y.  Mean  Value,  u  Standard  Deviation,  a 

1  Y1  Yi 


1-1 

0.121428e+04 

.1  u 

.  u 

1-2 

0.120000e+01 

.1  u 
.  Y2 

1-3 

0.358568e+02 

1  u 

.  y3 

1-4 

0.113389e+02 

.1  12 

1-5 

0.119523e-dl 

*2  UY5 

1-6 

0.10000Qe+01 

.2  u 

The  covariances  between  random  variable  pairs  in  y1>  1*1,6  are  given 
in  Table  4.2. 
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Table  4.2 

The  Covariance  Matrix 


.01.  u2x 

.001  ii£ 

0. 

0. 

0. 

0. 

i 

.01  u22 

0. 

0. 

°- 

0. 

•01  »23 

.001  ug 

.0001  jig  jig 

(symmetric) 

.01  »24 

0. 

0. 

i 

i 

' 

.04  U2g 

0  ... 

I 

■  i 

|  .04  „26 

The  aR/ari,  1*1,6  a<-e  computed  using  Equation  and  the  computer 


program  VARE.F.  The  results  are 

shown  in  Table  4.3. 

• 

* 

TABLE  4.3 

i 

Partial  Derivative  of  Peak  Response 

i 

3R/3Yj 

1  <k 

t 

a 

1*1 

0.57023e-02 

\ 

'  1*2 

-0.36344e+00 

i 

| 

1*3 

-0.17445e+00 

1*4 

-0.34095e-01 

■ 

1*5 

-0.13495e+02 

*< 

a 

1-6 

-0.42679e+0i 

i 

*  *. 

% 

Using  Equation  3-3  and  the  computer  program  VARE.F,  the  mean  of  maxi- 


man  displacement  is  computed.  The  results  Is: 

m 

•  ■»  2.65  (in)  I 

v 

Using  Equation  3-4  and  the  computer  program  VARE.F,  the  variance  of  % 

max: man  displacement  is  computed..  Each  term  in  Equation  3-4  is  listed  in  ^ 

Table  4.4.  j 

,  ,  C 

C 

v 

N 

9 
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Table  4.4. 

The  Terms  in  Variance  R 


.47943+00  |  -.3022-02  0. 


.1902e-02  |  0. 


.3913e+00  I  -.2413e-02 


,1494e~02 


(symmetric) 


.2017e-02 


,104ie-02 


.5339e-02 


0. 


,7268e+00 


The  variance  of  maximum  displacement  is  the  sum  of  every  term  as  shown  in 
Table  4.4.  The  result  is 

'/ar£ R]  *  1-617  (in^) 

Table  4.4  shows  that  only  the  variance  of  input  peak  air  pressure, 
natural  frequency  and  yield  displacement  have  significant  effect  on  the  • 
variance  of  maximum  displacement.  The  rest  of  the  terms  are.  almost 
negligible. 

4.2  The  Permanent  Offset  , 

In  Section  3.2,  the  permanent  offset  is  expressed  using  Equation  3-8. 
Let  H  represent  the  permanent  offset;  then 

H  »  h(3l,  3?*  33.  3<J 

inhere  H  «  permanent  offset, 

3^  ■  natural  frequency,  un 
32  *  VA7/m 

3^  ■  yield  displacement,  u7 
.  -if  *  maximun  displacement,  R 

The  means  and  standard  deviations  of  3^#  1*1,4  are  given  in  Table  4.5. 
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Table  4.5 

The  Means  and  Standard  Deviations  of  Si,  1*1,4 


Mean 


Standard  Deviation  a. 


1*1 

0.358569e+02 

-1  w8l 

1-2 

0.113389e+02 

**  ®S2 

1-3 

0.100000e+01 

•2  Ug3 

1-4 

0.265600e+01 

.12717 

Again  using  Equation  3-3  and  the  computer  program  VARE.F,  the  mean 
value  of  H  is  evaluated  and  the  result  is 

E[H]  *  1.4907  (in) 

The  first  order  derivatives  of  H  with  respect  to  bi,1*1,4  are  computed 
directly  frcn  Equation  (3-3)  and  ar*e  shown  in  Table  4.6. 

Table  4.6 

Partial  Derivative  of  H 

i  i-1  1*2  1*3  1*4 

3H/ 3g .  0.92387e-02  -0.29215e-01  -0.90000e+00  0.90000e+OC 


1-1 

0.92387e-02 


1-4 

0.90000e+00 


The  covariance  between  34  and  3^, 1*1,3  Is  computed  using  Equation 
(3-6).  The  covariance  between  two  3 j, 1*1,3  is  a  known  quantity.  There¬ 
fore,  when  the  variance  of  permanent  offset  is  expressed  by  the  form  of 
Equation  (3-4),  in  which  g  represents  permanent  offset  at  this  stage,  each 
term  in  Equation  (3-4)  is  computed  and  is  listed  in  Table  4.7. 

Table  4.7 

The  Terms  in  Variance  H 


.0011 


-  .00011 


.0011 


(symmetric) 


1.13100 


The  variance  of  H  is  the  sura  of  every  term  of  Table  4.7;  the  result  is 

Var  (H)  *  1.5906  (in2) 

From  Table  4.7,  it  is  observed  that  the  maximun  displacement  has  great 
influence  on  the  variance  of  permanent  offset.  The  covariance  between 
maximun  displacement  and  yield  displacement  has  less  influence,  and  the 
rest  of  the  terms  are  almost  negligible. 

4.3  The  Energy  Dissipated 

In  Section  3.3,  the  energy  dissipated  in  a  bilinear  hysteretic  system 
is  separated  into  two  parts.  The  first  part  is  energy  dissipated  in  the 
spring  and  the  second  part  is  energy  dissipated  in  the  damper.  Here,  only 
the  energy  dissipated  by  the  spring  is  considered.  From  Equation  3-9,  the 
energy  dissipated  by  the  spring  can  be  expressed  as 

£d  *  s  (ni»  n2»  na.  n4.) 
where  n1  *  n2  *  /A^/m,  03  *  u7,  ^4  *  R. 

The  means  and  variances  of  the  ni»  i*l*4  are  the  same  as  given  in 
Table  4.5.  The  first  order  derivatives  of  Ed  with  respect  to  ni,  i*l,4 
are  obtained  directly  from  the  expression  of  E^  which  is  given  in  Equation 
(3-10).  The  results  are 

Table  4.8 

The  Derivatives  of  Ed 

i  i*l  1-2  i-3  i-4 

3£d/ni  .8383728+03  -  .88 7303e+02  .56733e+03  .944165e+04 

The  mean  value  of  E^  is  computed  by  using  Equation  (3-3)  and  computer 
program  VARE.F  and  the  result  is  E(Ed)  *  .161578e+05  (lb-in).  When  using 
Equation  (3-4)  to  compute  the  variance  of  Ed,  the  variable  g  in  Equation 
(3-4)  represents  E^  at  this  stage.  Each  element  in  Equation  3-4  can  be 
obtained  using  the  approach  of  Section  4-2.  Table  4.9  lists  each  term  of 
Equation  (3-4). 


Table  4.9 

The  Terms  in  Variance  Ed 


.903687e+07 

-.302450e+05 

,34140Qe+G4 

-.131982e+08 

.101225e+05 

0 

.961467e+05 

'  .128976e+05 

-.921978e+06 

.144175e+09 

The  variance  of  Ed  is  the  sum  of  every  element  of  Table  4.9.  The 
result  is 

Var  (E^)  *  0.115133e+09  (lb2-  in2) 

From  Table  4.9,  it  is  observed  that  the  maximum  displacement  has  great 
influence  on  the  variance  of  E^,  and  the  rest  of  the  terms  are  almost 
negligible. 
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5.0  SUMMARY 

The  objectives  of  this  study  were,  first,  to  develop  a  model  that  . 
characterizes  the  permanent  offset  of  an  SDF,  bilinear  hysteretic  system 
subjected  to  blast  type  load;  and  second,  to  probabilistically  characterize 
the  features  of  tha  inelastic  structural  response.  In  this  report,  the 
response  characteristics  considered  were:  (1)  maxlmun  displacement 
response  (2)  permanent  offset  (3)  energy  dissipated  by  the  Inelastic 
spring.  The  statistical  properties  characterizing  these  measures  of 
response  are  the  mean  and  variance. 

In  Section  2,  the  maximum  displacement  response  was  computed.  The 
response  random  process  was  then  established.'  Section  3  discussed  the 
techniq-es  of  probabilistic  analysis  of  a  complicated  function.  Since  the 
response  random  process  which  was  established  in  Section  2  could  not  be 
expressed  in  a  closed  form,  the  computer  program  VARE.F  was  developed  to 
compute  the  moments  of  critical  measures  of  Inelastic  response. 

Several  numerical  results  were  shown  In  Section  4;  One  of  the  Impor¬ 
tant  results  showed  that  the  maxlmum^dlsplacement  response  has  great  Influ¬ 
ence  on  the  variance  of  permanent  offset  and  energy  dissipated.  The 
results  developed  In  this  study  are  restricted  to  blast  type  Inputs,  and 
bilinear  hysteretic,  SOF  systems. 

The  results  developed  in  this  report  can  be  used  In  the  probabilistic 
design  process.  For  example,  the  computer  program  VARE.F  may  be  used  to 
determine  the  mean  and  variance  of  some  critical  response  measures  and  the 
designer  may  determine  If  the  response  satisfies  certain  design  criteria. 
Or,  the  response  may  be  restricted  to  some  extreme  level,  then  the  system 
parameters  may  be  determined  by  using  computer  program  VARE.F  with  a  trial 
and  error  method.  Also  we  can  predict  the  mean  and  variance  of  the 
response  In  an  assessment  of  an  existing  building. 

Future  study  might  apply  the  present  techniques  to  multl-degree-of- 
freedom  systems.  Or  response  models  that  characterize  response  to  other 
than  blast  inputs  might  be  sought. 
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26 


VARE.F 


0*  *  *xx**'sr'jM»***ir5r**iiMr»r*'*iirx**iir*******ir*******»**x**ww*****x**w* 

C 

c ' 
c 

c  this  program  compute  the  var  [ED]  and  E  [ED] 
c 


c 

Equation  of  motion 

z 

c 

m7  y"  +  c7  y'  + 

k7 

(  )  »  fm  exp( -alpha  t  ) 

c 

or 

c 

y"  +  2  zita  wn 

y'  +  wn**2  (  y-0  )=A  exp  (-alpha  t) 

c 

where  : 

c 

g  is  z-max 

=G  (batal,bata2, . ,bata6) 

c 

b&tal=A, 

and  xmu(l)=E  [batal] 

c 

bata2=alpha 

$ 

and  xmu(2)=E  [bata2] 

c 

bata3=wn. 

and  xmu(3)=E  [bata3] 

c 

bata4=wnl. 

and  xmu(4)=E  [bata4] 

c 

bata5=zita. 

and  xmu(5)=E  [bata5] 

c 

bata6=v7. 

and  xmu(6)=E  [bataS] 

c 

c 

and  : 

c 

xrho ( i , j  ) 

are 

coefficient  of  correlation  between  bata(i) 

c 

xsigma( i ) 

are 

corresponding  STANDARD  DEVIATION 

c 

xvar  ( i , j  ) 

are 

vzriance  matrix 

c 

pgwb(i)  =d 

G  /d  bata(i) 

c 

c 

offset  is  the 

max 

offset  of  the  response, ie: 

c 

off set=H( gammal 

,  gamma2 , gamma3 , gamma4 ) 

c 

gamma 1=G 

and  zmul=E  [ gammal ] 

c 

gamma2=wyl 

and  zmu2=E  [ gamma2 ] 

c 

ganuna3=wn 

and  zmu3=E  [gamma3] 

c 

gamma4=v7 

and  zmu4=E  [gamma4] 

c 

and  zrho(i, j ) 

zsigma( i ) , zvar( i , j ) , ph( i )  represent  the  same 

c 

characteristic  as 

in  xmu 

c 

c  In  energy  dissipation  part  , 
c  smu(i)  are  related  to  the  spring 

c  therefore  :  , 

c  smu ( l)=E[wn] 

c  smu (2 )=E( wyl ] 

c  smu  ( 3  )  =E  [  v7  ] 

c  smu(4)=E[ zmax ] 

c  The  svar(i, j ) , srho(i, j ) ,pes( i )  are  corresponding 
c  variance  matrix'/  coeff .  of  correlation,  and  first  order 
c  ,  derivitive  w.r.p  to  those  parameter, 
c 

c  dmu(i)  are  related  to  the  damping 

c  Therefore  : 

c  dmu( 1 )=a 

c  dmu( 2 )=alpha 

c  dmu ( 3 ) =wn 

c  dmu(4)=zita 

c  dmu(5)=thita 


0  'kitieit'kicitirieitieiriririe'kitit'kicit'k’K’k'X'kTcifk  +  li'kit  Icic'xir-k'k-k'fcifkie'ririe-xieie-k'kieyc'k-x  x  x  *  *  ★  ★  x  * 
C 

dimension  xmu(6) ,xsigma(6) ,xrho(6,6) ,xvar(6, 6) ,pgwb(6) 
dimension  zmu(4) , zsigma( 4) , zrho(4, 4) , zvar(4, 4)  ,ph(4) 
dimension  smu(4) ,  ssigma(4)  ,srho(4, 4) ,  svar(4,  4)  ,pes(4) 
dimension  dmu(5) ,dsigma(5) 

dim. as ion  ff( 1024) . off set( 1024) , tryl( 1000} ,qql( 1000) ;v( 1024) 
real  ’*7,  m7 
c 

c  define  input  parameter 
fm=8500 . 
c7=6. 
m7=7 . 
k7=9C00 . 
a7*. I*k7 
v7=l . 
dt=0 . 01 
nb=1024 
a=fm/m7 
alpha* 1.2 
wn=(k7/ra7)**.5 
wyl=(a7/m7)**.5 
zita*c7/(2 .  *m7*wn) 
c 
c 

do  80  i=l,nb 
tx=(i-l)*dt 

ff  ( I  )*fm*exp(  -  alpha*  tx) 

30  continue 
c 

c  define  the  mbment  value  of  random  parameter 
xmu( 1 )=a 
xmu(2)«alpha 
xmu( 3 )=wn 
xmu( 4)=wyl 
xmu(5)=zita 
xmu ( 6 ) =v7 
c 

xsigma( 1)*. l*a 
xsigma(2)=0. l*xmu(2) 
xsigma(3 )*0. l*xmu(3 ) 
xsigma( 4)=0 . l*xmu( 4) 
xsigma<5)*.2*zita 
xsigma(6)*.2*v7 
c 

xrho(l, 2)=.l 

xrho <  3 , 4 )  = .  1 

xrho(3 ,  5  )*.  1 

xrho(3, 6)*. 01 

do  SO  1*1,6 

do  90  j*l, 6 

if  (i.eq.j)  go  to  91 


o  o  o  o  o  o 


xrho( j , i)=xrho(i,  j ) 
go  to  90 
91  xrho(i,j)=l. 

90  continue 
c 
c 
c 

c  find  the  value  of  q  which  is  the  parameter  to  be  identify 
c  in  offset 

call  sbilin  (FF, c7,m7, k7, a7, v7, dt, nb,v, offset, ened) 

thcon=of f set ( nb ) 

write  (6,10)  ened,  thcon 

10  format  ('Energy  dissioated  and  offset  from  sbilin  is  !, 
c  el3 . 6, 2x, f9 . 5 ) 


In  order  to  prevent  offset  oscillating,  check  if  the  max  offset 
equal  to  penitent  offset  ,  if  not,  it  means  that  the  forcinG' 
function  ,  or  the  system  itself  ,  does  not  practical. 

offmax=0. 
nbl=nb-l 
do  75  i*l,nbl 

offmax=amaxl  (offmax, offset (i ) ) 

75  continue 

if  (offmax. eq. 0. )  go  to  3000 
if  ( of finax.ne.  thcon)  go  to  3001 
.c 
c 

c  identify  the  parameter  ql  and  q2 
'  ql=1000. 
qx=ql 

o=thcon*2 . /3 . 1415926 
dif f=l . 
cycle=l . 
c 

1001  call  search  (offset, p,ql,nb,dt,epsiln) 
tryl( l)=epsiln 
ku=l 

cql(ku)=ql  , 

C 

dql=qx/( dif 5*10.7) 

C 

1021  qi=ql+dql 

call  search  (offset, p,ql,nb,dt,epsiln) 

ku=ku*l 

qql(ku)=ql 

tryl(ku)»epsiln 

if  ( tryl(ku) . It. tryl(ku-l) )  go  to  1021 
if  (ku.ge.3)  go  to  2000 
c 

ql=ql-dql 
1C  10.  ql-ql-dql 


:t/P,ql,nb,  dt,  epsiln) 
yl(ku-l) )  go  to  2000 

JO  TO  1001 


^0}********************************** 

^l^rvt******************************* 

•*7thlr)  3  ' *WU* 4) ' x:au( 5 > ' «“(6). 

,i )  ,xslgaa( i ),!«!, 6) 


!  58  zmax  is  \fl6.9) 

*3  :  , 6x, ' X3igna( i J  is  ') 


it, k7, a7,dt,nb,pmox,pg) 


3(  i ) ,  i*l,  5) 
»ta  is  ;  ' ) 


xsigmafiqq) 

xrho(ipp,iqq) .eq.O. ) 
3  to  98 
3  to  97 


go  t6  99 


o  o 


pga=pgwb ( i cn ) 

icn=iqq 

pgb=pgvb(icn) 

xvar( ipp, iqq)=pga*pgb*sigx*xrho( ipp, iqq) 
go  to  100 

97  icn=ipp 
pga=pgwb( icn) 

xvar( ipp, iqq)=(pga) **2*sigx*xrho( ipp, iqq) 
go  to  100 

98  xvar(ipp, iqq)=xvar(iqq, ipp) 

go  to  100 

99  xvar(ipp, iqq)=G. 

100  continue 


do  200  i=l/6 
do  200  j=l,6 
varx=varx+xvar ( i , j ) 

200  continue 

write  (6,24) 

write  (6,25)  ( (xvar (i, j ) , j=l, 6) , i-1, 6) 
write  (6,26)  varx 

24  format  (/' X-variance  matrix  :  ') 

25  format  (5el6.6) 

26  format  (/'  var  [G]  is  ' ,2x,ei3.6) 
c 

c 

c 

c 

c 

c  compute  E  [offset! 

off=(pmr-v7 ) * ( 1 . -a7/k7) 
c 

c  .  ,  . 

£  ********************************************  ***************** 
C 

c 

c 

c  COMPUTE  VAR  [offset] 
c 

c  ph( l)»dh/d( gamma 1) 
c  ph(2)»dh/d(gamma2 ) 
c  ph( 3 )=dh/d( gamma3 ) 
c  ph( 4)=dh/d(gamma4) 
c 

c  ************************************************************* 

c 

c 

ph( 1)»1 . -(wyl/wn)**2 
ph(2)*(pmr-v7)*(-2. )*wyl/wn**2 


c 


ph(3 )=(pmr-v7 ) *2 . *wyl**2/wn**3 
ph(4)=-l . * ( 1 . - (wyl/wn) **2) 


zmu ( 1 ) =pmr 
zmu(2)=wyl 
zmu ( 3 ) =wn 
zmu(4)=v7 


zsigma( 1 )=varx** . 5 
zsigma(2)=xsigma( 4) 
zsigma( 3 )=xsigma( 3 ) 
zsigma( 4)=xsigma( 6) 
write  (6,55) 

write  (6,56) (zmu( i ) , zsigma( i ) , i-1,4) 
write  (6,27)  off 

27  format  (/'S  (offsetl  *  ',fi5.3) 

55  format  ('  zmu(i)  is  : ' , 6x, ' zsigma  is 

56  format  (el4. 6, 3x, el4. 6) 


c  define  zrho(i,j) 

zrho (2,3) =xrho (3,4) 
zrho (3,4) =xrho (3,6) 
do  320  i»l, 4 
do  320  j=l, 4 
if  (i.eq. j )  go  to  321 
zrho( j , i )=zrho( i, j ) 
go  to  320 
321  zrho( i, j )»1. 

320  continue 


c 

c  1 

c  cbrapute  covfG  wyl I ,  cov(G  wn],- covfG  v7],  E(G**2] 
1  do  350  i*l,4 

if  (i.eq. 1)  go  to  351 

call  egr  (xmu, xsigma, xrho, pgvo, smu ( i ) , epara ) 
zvar( 1, i )*epara#ph( 1) *ph( i) 
go  to  350 
351  zvar(l,i)*(zsigma(l)*ph(.l)  )**2 
350  continue 


c 

c 


do  400  1*2 , 4 
do  400  j*2,4 
if  (i.gt.j)  go  to  400 

zvar( i, j )*zrho( i , j )*zsigma( i )*zsigm»( j )*ph( i)*ph( j ) 
400  continue 

do  410  1*1,4 
do  410  j*i ,  4 
zvar( j , i )*zvar( i, j ) 

410  continue  . 
write  (6,28) 


write  (6,29)  (ph( i ) , i=l, 4) 
write  (6,30) 

write  (6,31)  ( (zvar(i, j ) , j=l, 4) , i=l, 4) 

28  format  (/'ph(i)  is  :  '  ) 

29  format  (el6.7) 

30  format  (/'2-variance  matrix  is  :  '  ) 

31  format  (4fl5.6) 
c 

varof=0.  • 

do  500  i=l,4 
do  S00  j=l ,  4 
varof=varof+zvar( i , j ) 

500  continue 
c 

write  (6,32)  varof 

32  format  (/'Var  [offset]  =  ',fl5.3) 
c 

c 

w  ★  ★  #  *  it  it  ★  *  *  »»  ***#**#*'****★**★***★★★★  w*  ★»*★%*★**★****★★  ★ 

C 

c 

c 

c  compute  E  [ED] 
c  ED  :  ENERGY  Dissipated 

c  ED1=  dissipate  by  spring 

c  ED2=  dissipate  by  damping 

c 

£'«r*****»****'#Tir'**'****i»»********#**i«**llMHr************************ilr 

C 

c 

c  define  raeam  value  and  standard  deviation 
smu ( 1 ) =vn 
smu(2)=wyl 
smu ( 3 ) -v7 
3ir.u(  4)=pmr 
c 

ssigma(  1  )=>xsigraa(  3 ) 
ssigma(2 )=xsigma( 4) 
ssigma(3 )=xsigma( 5) 
ssigma(4)*varx** . 5 
c 

dmu(l)=a 
dmu(  2 )  =*aipha 
dmu ( 3 ) =wn 
dmu(4)=zita 
dmu(5)=off 
c 

dsigma( 1 )=xsigma( I ) 
dsigma(2 )=xsigma(2 ) 
dsigma(3 )axsigma( 3 ) 
dsigma(  4)s*xsigraa(  5 ) 
dsigma( 5 )=varof . 5 


C3.ll  engl ( smu, m7, edl ) 
call  eng2  (dmu,  nb,  q,  m7,  ed2  ) 
write  (6,50) 

write  (6,61)  ( snu( i ) , dmu( i ) , i=l, 4) 
wri te  (6,62)  dmu ( 5 ) 
write  (6,33)  edl,ed2 

60  format  ('  smu(i)  is  ',8x,'dmu(i)  is  ') 

61  format  (2el5.6) 

62  format  (15x,el5.6) 

33  format  (/'ElEdlJ  is  ',615.6,'  E [ Ed2  ]  If  ’,el3.6) 

/*• 

c 

c 

c 

c 

c  compute  var  (Ed] 
do  600  cn=l,4 
call  pengl(cn, smu,m7,pedl) 
icn=cn 

pes( icn)=pedl 


600 

continue 

write 

(6,34) 

write 

(6,35)  (pes( i ) , i=i, 4) 

34 

format 

(/'oes(i)  is  : ' ) 

35 

format 

( e!3 . 5 ) 

c  define  3rho( i , j ) , j=l , 3  1=1,3 

c  where  srho(i,j)  is  the  coefficient  of  correlation 
c  in  spring  snergy  dissipate  case 
c 

srho ( 1 , 2 ) =xrho ( 3 , 4 ) 
srho ( 1 , 3 ) =xrho ( 3 , 5 ) 
do  650  i=l, 4 
do  650  j=l ,  4 
if  (i.eq.j)  go  to  651 
srho( j , i )=srho( i , j ) 
go  to  650 
651  srho(i,J)=l. 

650  continue 


do  700  i=l, 3 
do  700  j*l, 3 
if  ( i . gt. j )  go  to  700 

svar(  i ,  j  )=srho(  i,  j )  ,l,ssigna(  i  )*3sigma(  j  )  *pes(i )  *pes(  j  ) 
continue 

svar( 4, 4)»/arx*pes( 4)**2 
do  300  i*l, 3 

call  egr(xmu,  xsigma,  xrho,  pg*vb,  smu(  1) ,  covhs) 
svar( i , 4) =covhs#pes( i ) *pes( 4) 


7C0 

c 

c 


34 


o  n 


800  continue 


do  900  1=1,5 
do  900  j=l, 5 
svar( j , i)=svar(i, j ) 

900  continue 

write  (6,36) 

write  (6,37) ( ( svar( i, j ) , j-l, 4) , i=l , 4) 

36  format  (/' S-variance  matrix  is 

37  format  (4el5.6) 
c 

do  920  i=l, 4 
do  920  j*l,4 
svarl=svar 1 +svar ( i , j ) 

920  continue 


write  (6,38)  svarl 
38  format  (/'var  [EDI)  IS  ’,el3.6) 
go  to  5000 

3000  WRITE  (6,18) 

18  FORMAT  (’  MO  YIELDING  OCCUR’) 

GO  TO  6000 

3001  WRITE  (6,19) 

19  FORMAT  (  ’  THIS  CASE  DOEN  NOT  PRACTICAL’ ) 

6000  STOP 
END 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


Thd 


subroutine  <pgl( icn, xmu, ff , k7, a7, dt,nb,pmo,pg; 
s  subroutine  .ompute  dG/db,  denote.it  as  pg  ; 
’random  parameter. 

1  number  ; 


Where  b  is 
cn  =  contrc) 
if 


cn=l. 

=2. 

=3. 

=4. 

=5. 

=6, 


compute  dG/dA 
conipute  dG/dalpha 
compute  dG/dwn* 
compute  dG/dwyl 
conpute  dG/dzita 
conpute  dG/dv7 


dimension  q 
real  K7,M7, 

a=xmu ( 1 ) 
alpha=xmu(^ 
wn=xmu ( 3 ) 
wyl=xmu( 4) 
zita=xmu( 5 ) 


f (1024) ,xmu(6) 
k71,  Jc72 


) 


"j-V 


< 


v7=xmu ( 6 ) 

rc  ad  (5, * )  error 
error=500. 

if  (icn. eq. 6) , go  to  6 
if  (icn. eq. 5).  go  to  5 

if  (icn.eq.4)  go  to  4 

if  (icn.eq.3)  go  to  3 

if  (icn.eq.2)  go  to  2 

if ' (icn. eq. 1 )  go  to  1 
go  to  4999 

m7=k7/(wn**2 ) 

fm=a*m7 

dfm-f m/error 

da=dfm/m7 

al=a+da 

a2=a-da 

call  g( al, alpha, wn, wyl, zita, v7,  ff, k7, a7,dt,nb, spmo, thir) 
write  (6,22)  thir 

gl-spmo 

call  g(a2 , alpha, wn, wyl, zita, v7 , ff , k7, a7, dt,nb, spmo, thir ) 
g2=spmo 

pg=(gl-g2)/(2.*da) 

pmo=gl 

go  to  1000 

dalpha=alpha/error 

alphal=alpha+dalpha 

alpha2=alpha-dalpha 

call  g( a,  alphal, wn, wyl, zita, v7, ff , k7, a7, dt, nb, spmo, thir) 
write  (6,22)  thir 

gl=spmo 

call  g( a, alpha2 , wn, wyl, zita, v7, ff , k7,a7,dt, nb, spmo, thir) 
g2=3pmo  , 

pg=( gl-g2 )/ ( 2 . *dalpha ) 
pmo=gi . 
go  to  5000 

ak7=k7/error 

k71=k7+dk7 

k72=k7-dk7 

m7=a7/(wyl**2 ) 

wnl=(k71/m7)**.5 

wn2=(k72/ra7)** . 5 

dwn=wnl-wn2 

call  g( a, alpha, wnl, wyl, zita, v7, ff , k71, a7, dt, nb, spmo, thir) 
write  (6,22)  thir 

gl=sprao 


on  f>'  o  o  rfi  o 


call  g( a, alpha, wn2, wyl, zita, v7, ff , k72, a7,dt, nb, spmo, thir) 
g2=spmo 

pg=(gl-g2)/(2.*d’-m) 
pmo=gl 
go  to  5000 

da7=a7/error 
a71=a7+da7 
a72=a7-da7 
m7=k7/(Wn**2 ) 
wyll=(a7l/m7)**.5 
wyl2=(a72/m7)**.5 
dwyl=wyll-wyl2 

call  g( a, alpha, wn,wyll, zita, v7, ff, k7, a71,dt, nb, spmo, thir ) 
write  (6,22)  thir 

gl=spnu> 

call  g( a, alpha, wn, wy!2 , zita, v7, ff , k7, a72, dt, nb, spmo, thir ) 
g2=spmo 

pg=(gl-g2)/(2.*dwyl) 
pmo=gl 
go  to  5000 

5  dzita=zita/error 
zital»zita+dzita 
zita2=zita-dzita 

call  g(a, alpha, wn, wyl, zital, v7, ff , k7, a7, dt,nb, spmo, thir ) 
write  (6,22)  thir 

gl=spmo 

call  g(a, alpha, wn, wyl, zita2, v7, ff , k7, a7,dt,nb, spmo, thir) 
g2=spmo 

pg=(gl-g2)/(2.*dzita) 
pmo=gl 
go  to  5000 

6  dv7=v7/error 
v71=v7+dv7 
v72=v7-dv7 

call  g(a,  alpha, wn, wyl, Zita, v71,ff,k7,a7.,dt,nb,  spmo,  thir) 
c  write  (6,22)  thir 

c 

gl=spmo 

call  g(a, alpha, wn, wyl, zita, v72, ff ,k7, a7,dt, nb , spmo, thir) 
g2-spmo 

pg=(gl-g2)/(2.*dv7) 
pmo=gl 
go  to  5000 
c 

22  format  ( ’ thita  ratio  is  ’ , fl3.7) 

4999  write  (6,19) 

19  format  ('Control  Number  Wrong') 
c 

1  \ 

\. 

\ 


5000  return 
end 
c 
c 

subroutine  g( a, alpha, wn, wyl , zita, v7, ff , k7, a7, dt, nb, pmr, thir) 
C  This  subroutine  is  to  predict  the  max  offset  for  blast  load 
C  pmo=predicted  max  response 

c  pmr=100:  represent  jmm  expired  (in  do  loop) 
t  pmr=101:  represent  something  trouble  in  finding  tmax 
c 

dimension  ff(  1024) ,  W(  1024) ,  0FFSET(  1024) ,  SLOPE(30) 

REAL  K7,M7 

WD-( ( 1 . -Z1TA**2 ) ** . 5 ) *WN 
m7=k7/'(wn**2 ) 
c7=2 . *zita*wn*m7 
V=-ZITA*WN 
F=ALPHA+V 
C 
C 

wy=(wyl**2-(zita*wn)**2)** .5 
D=(WN**2-WY1**2 ) *V7 
C 
•  c 

call  sbilin  ( ff,  c7,  m7,  k7,  a7,  v7,  dt,  nb,  w,  offset,  ened) 
c 
C 

C  CALCULATE  T7 

B=2 . *AL?HA/( WD**2+AL?HA**2 ) 

C=2 . *77/ A 

T7=( B+ (3**2+4.*C)**.5)/2. 

T71=T7 

c  WRITE  (5,1)  T7 

c  1  FORMAT  (’THE  TEMPERATORY  T7  =  ’  ,  Fll.  7) 
c 

C  CALCULATE  EXACT  T7 
•  DO  700  JI 1*1, 2000 
XX=JII*.00001 
T7=T71+XX 

■  ■  Xl=(  (ALPEA-ZITA*WN)**2.+WD**2)*WD 
X2=A/X1 

X3=EXP  ( -ZITA*WN*T7 ) 

X4=(ALPHA-ZITA*WN)*SIN'WD*T7) 

X5=WD*C0S(WD*T7) 

X6=WD*EX? ( -AL?HA*T7 ) 

C 

C  CHECK  Z  AT  T7  IF  IT  IS  EQUAL  TO  77 
ZT7=X2*  ( X3*  (X4-X5  )  +X6 ) 

ERR1=A3S ( 1 . -ZT7/V7 ) 

IF  (ERR1.LE. 0.0005)  CO  TO  701 

700  CONTINUE . 
c 

c  compute  the  z 1 ( t7 ) 

701  X7=WD  * ( ALPHA- ZITA* WN ) * COS ( WD *T7 ) 


X8=WD**2*SIN(WD*T7) 

X9=-ZITA*WN*X3* (X4-X5 ) 

ZD0T7=X2* (X9+X3* (X7+X8 ) -ALPHA*X6 ) 

C 

*  c 

c  calculate  tmax  (using  approached  method) 
rx=2 . *zita*wn 
gl=a* ( 1 . -exp ( - rx*  t7 ) ) 
g2=(wn**2)*v7*(exp(-rx*t7)-l. ) 
gx=- ( gl+g2 ) +zdot7*rx 
c 

hl=a*exp(alpha*t7) 
h2=a*rx/( alpha- rx)+(wn**2 ) *v7 
h3 =a * alpha* exp ( t7* ( alpha- rx) )/( alpha- rx) 
h4=gx*exp(rx*t7) 
h5=hl-h3 
h6=h2+h4  ' 
c 

pl=h5*alpha**2+h6*(rx**2) 
p2  =h5  *  a lpha +h6  *  rx 
p3=h5+h6-(wn**2)*v7 
c 

ux=p2**2-4. *pl*p3 
.  ,  if  (,ux.  le.O. )  go  to  2700 

uu=ux** . 5 

tmaxl=(p2+uu)/(2 . *pl ) 
tmax2=(p2-uu)/(2 . *pl) 

if  ( tmaxl . le.O. 17. and. tmaxl. ge . 0. )  go  to  70 

tmax=tmax2 

go  to  71 

70  tmax=tmaxl 

71  tmaxt=tmax 
c 

c  1 

c,  compute  cl  and  c2  . 

ql=-d/(wyl**2 ) 
c 

qxl=alpha**2 
qx2=-2 . *zita*wn* alpha 
qx3=wyl**2 
q2=a/(qx4.+qx2*qx3 ) 
c 

xal=exp( -zita*wn*t7 ) 
xa2=-zita*wn*cos(wy*t7) 
xa3=s-wy*3in(wy*t7) 
xat7=xal* (xa2+xa3 ) 
c 

xbl=-zita*wn*sin(wy*t7) 
xb2=wy*cos(wy*t7) 
xbt7=xal* (xbl+xb2 ) 


xct7=*-alpha*q?*exp(  -alpha*t7) 


c 


no  o  o 


xdt7=xal*cos ( wy*t7 ) 

xet7=xal*sin(wy*t7) 

xft7=ql+q2*exp(-alpha*t7) 

xgt7=v7-xft7 

xht7=zdot7-xct7 

deno=xat7*xet7-xdt7*xbt7 
cl=(xht7*xet7-xbt7*xgt7)/deno 
c2=( xat7*xgt7-xdt7*xht7 )/deno 

compute  exact  tnax,  since  z'=0  when  t=tmax 
tmx=tmaxt 
dtx-tmaxt 

880  tx=dtx/30. 

do  900  j mm=l , 30 
tmx=tmx-tx 

xal=exp( -zita*wn*tmx) 
xa2=-zita*wn*cos(wy*tmx) 
xa3 =- wy *  s i n ( wy *  tmx ) 
xa ttax=xal * ( xa2  +xa3 ) 


xbl=-zita*wn*sin(wy*tmx) 
xb2 =wy * c o s ( wy * tmx ) 
xbtmx=xal* (xbl+xb2 ) 

xctmx=*-alpha*q2*exp( -alpha*tmx)  1 

z d tmx=x a tmx * c 1 + xb tmx * c 2 + x c t mx 
slope( jmm)-zdtmx 

if  ( jmm. eq. 1 . and. slope( jmm) . gt . 0. )  go  to  899 
if  (jmm.eq.l)  go  to  900 
proa=slope( jmm) *slope( jmm-1 ) 
if  (prod.le.0.)  go  to  902 
900  continue 

write  (6,10) 

10  format  ('something  WRONG  in  t-m4x( 3ubG] ' ) 

9C2  if  ( slope( jmm-1 ). gt. -. 0001 )'  go  to  901 
dtx-dtx/30. 
tmx=tmx+tx 
go  to  880 

899  do  920  jmm=l,1000 
tx=( jmm-1) *0.0001 
tmx=tmaxt+tx 
xal=exp(-zita*wn*tmx) 
xa2=-zita*wn*cos(wy*tmx) 
xa3  *- wy*  sin ( wy *  tmx ) 


xatmx=xal*(xa2+xa3 ) 
c 
c 

xbl=-zita*wn*sin(wy*tmx) 
xb2=wy*cos(wy*tmx) . 
xbtmx-xal* ( xbl+xb2 ) 
c 

xctmx--alpha*q2*exp( -alpha*tmx) 
c 

z d tmx=x a tmx * c 1 + xb tmx * c 2 + x c tmx 
if  (zdtmx. le.0.01)  go  to  901 
920  continue 

go  to  2100 
c 

c  compute  exact  z(at  exact  tmax) 

901  xdtmx=xal*cos(wy*tmx) 
xetmx=xal*sin( wy*tmx ) 
xftmx=ql+q2*exp( -alpha*tmx) 
ztmx=xdtmx*c l+xermx*c2 +xf tmx 
c 

c  calculate  predicted  max-offset 
pmo=( ztmx-v7 ) * ( 1 . -a7/k7 ) 
thi r=pmo/o f f se t ( nb ) 
pmr=ztmx 

go  to  3000 

* 

write  (6,98.) 
format  ( ' fmm  expire* ) 
pmr=100. 
go  to  3000 
write  (6,99) 

format  ('U  less  than  zero') 
pmr=101. 
return 
END 

SUBROUTINE  S3ILIN( F, C7 ; M7 , K7 , A7 , 77 , D9 , N , V, V0 , aned) 
DIMENSION  F(1024) , V0 ( 1024) ,V( 1024) 

DIMENSION  VI ( 1024), V2 ( 1024) 

REAL  M7,K7,K8,K9 
U7=K7*V7 
K9=l . -A7/K7 

C  INITIALIZE  VARIABLES 
V(1)»0. 

V0(1)=0. 

V1(1)=0. 

V2(1)*F(1)/M7 

C  START  THE  RESPONSE  CYCLE 
Q1=6.*M7/D9**2 
Q2=3.*C7/D9 
Q3=6 . *M7/D9 


c 

c 

2100 

98 

2700 

99 

3000 

C 


:< 


i 


8 


* 


■ 


Q4=3 . *M7 
Q5=3 . *C7 
Q6=0 . 5*D9*C7 
Q7=3./D9 
Q8=3 .  ■ 

Q9=.5*D9 

K8=K7 

NM=N-1 

ENED-O. 

DO  1199  1=1, NM 
11=1+1 
U1=Q1+Q2+K8 
U2=Q3  * VI ( I ) +Q4* V2 ( I ) 

U3=Q5*V1 ( I ) +Q6*V2 ( I ) 

V5=(F( 1 1 ) —  F ( I )+U2+U3 )/Ul 
V6=Q7*V5-Q8*V1(I)-Q9*V2(I) 

V( II )=V( I ) +V5 
V1(:I1)=V1(I)+V6 

C  COMPUTE  THE  STIFFNESS  AT  T+DT 
X0=K7* ( V ( 1 1 ) - VO ( I ) ) 

X1=A7* ( V ( I 1 ) - V7 ) +U7 
X2=A7  * ( V ( 1 1 ) +V7 ) -U7 
IF  (XO.GT.X1)  GO  TO  1150 
IF  (X0.LT.X2)  GO  TO  1160 
K8=K7 

V0(I1)=V0(I) 

GO  TO  1170 

1150  IF  (Vl(U)  .GT.O.  )  K8=A7 
IF  (Vl(Il).LS.O.)  K8=K7 
V0(I1)=(V(I1)-V7)*K9 
GO  TO  1170 

1160  IF  (71(11) .LT.O. )  K8=A7 
IF  (Vl(Il).GS.O. )  K8=K7 
V0(I1)=(V(I1) +V7) *K9 

1170  V2(I1)  =  (F(I1) -C7*V1( I1)-K7*(V(I1) -V0( II) )  )/M7 
ened=ened+d9* i 5  * ( vl ( i ) * ( f ( i ) -m7  * v2 ( i ) )  + 
c  vl(il)*(f(il)-m7*v2(il))) 

1199  CONTINUE 
RETURN 
END 
c 
c 

subroutine  fft  { *, n, nb, isgn, dt ) 
complex  a(nb),u,w,t 
c 

c.  dividing  all  element  by  nb 

do  1  j=l,nb 
1  ag)=a(j)/nb 
c 

c  reorder  sequencr 

nbd2=nb/2 
nbml=nb-l 
J=1 


42 


*r 


do  4  1=1 , nbml 

<  if  (l.ge. j )  go  to  2 

t=a(j) 
a(j)=a(l) 

I  a ( 1 ) =t 

2  k=nbd2 

3  if  (k.ge.j)  go  to  4 

j=j-k  . 
k=k/2 
go  to  3 

4  j-j+k 

c  calculate  fft 

pi=3. 1415926 
do  6  m=l , n 
u=(l. 0,0.0) 

ME=2**M 

k=me/2 

w=cniplx(cos(pi/k)  ,  isgn*sin(pi/k) ) 

do  6  j=l,k 

do  5  l=j,nb,me 

lpk=l+k 

t=a( lpk) *u 

a( lpk)=a( l)-t 

5  a(L)=a( l)+t 

6  u=u*w 
c 

tt=(nb-l ) *dt 

»  c 

if  (isgn.eq.l)  go  to  99 
do  110  i=l, nb 
a(i)=a(i)*tt 
110  continue 
go  to  100 
c 

99  do  120  i*l,nb 
a(i)=a(i)/dt 

120  continue 

100  return 
end 

c  _ 

c 

c 

subroutine  engl  ( smu , ra7 , endl ) 
c  this  subroutine  compute  the  S  [EDI] 
dimension  smu(4) 
real  k7,m7 
k7=( smu( 1)**2 )*m7 
a7=(smu(2)**2)*m7 
v7=smu(3) 
zmax=smu(4) 

thita=( 1. -a7/k7)*(zmax-v7) 
c 

xl*. 5*k7*v7**2 


x2=a7* ( zmax-v7 ) **2/2 . 

<  x3=k7*v7*(znax-v7) 

x4=k7* ( zmax-thita) **2/2 . 
c 

%  endl=xl+x2+x3+x4 

return 
end 
c 
c 
c 

subroutine  eng2  ( dmu,  nb,  q,m7,  ed2  ) 
c  this  subroutine  compute  E  [ED2] 

dimension  dmu( 5 ) , fl ( 1024) , £2 ( 1024) , w( 1024) , deno(512 ) 
complex  ffl( 1024), ££2( 1024), ff( 512 ),hh(512) 
real  ra7 
c 

a=dmu ( 1 ) 
alpha^dmu ( 2 ) 

wn=dmu ( 3 )  . 

zita=dmu(4) 
thita=dmu(5) 
p=thita*2 . /3 . 1415926 
c 

dt=0 . 01 
tt=(nb-l)*dt 
c 

do  100  i»l, 1024 
*  tx=( i-1 ) *dt 

fl(i)=a*exp(-alpha*tx) 

£2(  i)s=wn**2*p*atan(q*tx**4) 
ffl( i )=cmplx  (fl(i) ,0. ) 
ff2(i)=cmplx  ( f 2 ( i ) , 0 . ) 
w( i )=( i-1 ) *2 . *3 . 1415925/tt 
100  continue 
c 

call  f£t(f£l, 10,nb,-l,dt) 
call  fft( f£2, 10,  nb,  -l,d,t ) 
c 

nb2=nb/2 

do  200  i»l,nb2 

ff(i)=f£l(i)+ff2(i) 

re*wn**2-w(  i,)  **2 

xia*-2 . *zita*wn*w( i ) 

deno( i )»re**2^xim**2 

hh(  i  )=*cmplx  ( re, xim)/deno( i ) 

200  continue 
c 

end2=0 . 

do  400  i»l, nb2 
xel*cabs(hh(  i ) ) 
xe2*cabs( ff (i ) ) 
ed2*ed2+(xel*xe2*v( i ) )**2 
400  continue 


2 

ff 


e<±2=2 .  *ed2 
return 
end 
c 
c 
c 

subroutine  pengl(cn, smu,m7,pedl) 
c  this  subroutine  compute  d  EDI/d  bat a 
c  cn  =1  d  EDI/d  wn 

c  cn  =2  d  EDI/d  wyl 

c  cn  =3  d  EDl/dv7 

c  cn  -4  d  ED 1/d  zmax 

c 

dimension  smu(4) 
real  m7 
wn=smu( 1 ) 
wyl=smu(2) 
v7=smu(3) 
zmax=smu(4) 

const=( zmax-v7 ) *wyl**2/wn+v7*wn, 
c 

if  (cn.eq.l.)  go  to  100 

if  (cn.eq.2.)  go  to  200 

if  (cn.eq.3.)  go  to  300 

if  (cn.eq.4.)  go  to  400 

go  to  600 
c 

ICO  constl=const*  ( v7-  ( zmax-v7 )  *  ( wyl/wn )  **2 ) 

pedl=( wn*v7**2*2 . *vn*v7* ( zmax-v7 ) -constl ) *m7 
go  to  700 

200  pedi=( wyl* ( zmax-v7 ) **2-2 . * ( zmax-v7 ) *wyl*const/wn) *m7 
go  to  700 

300  const2=const* ( wn**2-wyl**2 )/vn 

pedl=( zmax-v7 ) * ( wn**2-wyl**2 ) -const2 
go  to  700 

400  pedl=( (zmax-v7)*wyl**2+v7*wn**2-const*wyl**2/wn)*ra7 
go  to  700 

600  write  (6,1) 

1  format  ('  control .number  wrong’) 

700  return 
end 
c 
c 
c 
c 

c  . 

C  ' 

subroutine  egr  (xmu,xsigna,xrho,pgwb, gamma, epara) 
c  this  subroutine  compute  E[R, gamma] 
c  where  R=zmax 
c  xau:  mean  value  of  bata 

c  where  bata  see  main  program 


c  pgwb:  d  g/d  bata(i) 
c 

dimension  xmu(6) ,xsigma(6) , xrho(6, 6) ,pgwb(6) 
c 

do  100  i=l , 6 

if  (xmu( i ). eq. gamma)  go  to  101 

100  continue 
write  (6,1) 

1  format  ('gamma  is  not  equal  to  xmu(i)  in  subroutine  egr, WRONG 

101  nj=i 
c 

epara=0 . 
do  200  3=1,6 

exl=xrho(nj , j )*xsigma(nj )*xsigma( j ) 
ex2=(exl)*pgvb( j ) 
epara=epara+ex2 
200  continue 
return 
end 
c 
c 
c 

subroutine  search  ( thita, p, a, nb, dt , epsiln) 
dimension  thita(rJb) 
sum 1=0 . 
do  100  i=l,nb 
tx=(i-l)*dt 
ty=p*atan(a*(tx**4) ) 
suml=suml+ ( thita( i )- ty ) **2 
100  continue 

epsiln=suml 

return  , 

end  .  | 


